A Solution Routine for Singular Boundary Value Problems

In this report, we discuss the implementation and numerical aspects of the Matlab solver sbvp designed for the solution of two-point boundary value problems, which may include a singularity of the first kind, z′(t) = f(t, z(t)) := 1 (t− a) · z(t) + g(t, z(t)) , t ∈ (a, b), R(z(a), z(b)) = 0. The code is based on collocation at either equidistant or Gaussian collocation points. For singular problems, the choice of equidistant nodes does not imply a loss of efficiency since no convergence order higher than the stage order can be expected in general. An error estimate for the global error of the approximate solution is also provided. This estimate is obtained by a modification of the Defect Correction idea originally proposed by Zadunaisky in 1976. This estimate has been proven to be asymptotically correct and provides the basis for an adaptive mesh selection strategy. Here the grid is modified with the aim to equidistribute the global error. Most importantly, we observe that the grid is refined in a way reflecting merely the smoothness of the solution, unaffected by the singularity of f . We discuss details of the efficient implementation in Matlab 6 and illustrate the performance of the code by comparing it with the standard Matlab solver bvp4c and the Fortran 90 code COLNEW. The sbvp package is available from http://www.math.tuwien.ac.at/~ewa.